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We present an analytic calculation of the gravitational field inside and near a caustic ring of dark 
matter. The calculation may facilitate N-body simulation studies on the effects that dark matter 
caustics have on galaxy formation. 

I. INTRODUCTION 

One mystery in modern cosmology is the identity of the dark matter, which constitutes 23% of the energy density 
of the current universe. By now, a host of observations as well as numerical simulations of structure formation 
have established that dark matter must be cold and collisionless nonbaryonic particles that interact primarily via 
the gravitational interaction. Since no known particles have these properties, one has to turn to physics beyond the 
Standard Model for suitable candidates. Currently, the leading candidates of cold dark matter (CDM) are axions and 
weakly-interacting massive particles (WIMPs). The former originates from Peccei and Quinn's dynamical solution to 
the strong CP Problem [IH1], while the latter naturally arises from R-parity-conserving supersymmetry. 

For a long time, it is believed that the two candidates are indistinguishable on observational grounds, even though 
fundamentally they are very different. The reason for this has to do with the fact they can both be viewed as 
classical approximations of quantum fields, albeit in different limits. WIMPs, as fermions, arc in the regime of 
classical particles, in which we take h — > and to, k — > oo while keeping their energy E — huj and momentum p = hk 
finite. On the other hand, axions, as bosons, reside in the classical field regime, in which we take h — > and the 
occupation number N — > oo for constant ui and k to keep E = Afhuj and p — Nhk fixed [5] . Hence, both behave as a 
pressureless perfect fluid, and have similar properties as far as large scale structure is concerned. For example, it has 
been demonstrated that the evolution of cosmological perturbations in an axion BEC and ordinary cold dark matter 
are practically indistinguishable on scales of observational interest provided the BEC does not rethermalize [BHH] • As 
such, it is widely held that the verdict on the true identity of dark matter can come only with a direct detection 
in experiments (there is of course the possibility that neither candidate is correct). Currently, a large number of 
experiments worldwide are actively searching for both candidates, so far with negative results. 

Recently, however, it has been pointed out that dark matter axions, being a Bose-Einstein condensate (BEC) [7J, 
could potentially give rise to novel avenues to differentiate them from WIMPs. If, for example, thermal contact 
were established between the axions and photons after big-bang nucleosynthesis and before recombination, various 
cosmological parameters (namely the effective number of neutrinos N e ff and the primordial abundance of light 
elements such as 4 He, 3 He, 7 Li, etc.) would change [9]. Hence, a more precise measurement of N e ff by the Planck 
mission in the near future might help shed light on the identity of dark matter. 

In addition to cosmological parameters, the study of the catastrophe structure of inner caustics of galactic halos may 
provide yet further clues |10j . Caustics are formed as galaxies accrete dark matter particles in their vicinity during 
structure formation. They are surfaces of high density which can be viewed as the envelope of particle trajectories, 
as dark matter particles fall in and out of the gravitational well. This process of particle infalling and outfalling 
occurs repeatedly, and a set of inner (as well as outer) caustics form. It turns out that there exists a dichotomy in the 
classification of the inner caustics in terms of their catastrophe structure, depending on the distribution of angular 
momentum of the infalling particles. The first case is characterized by "net overall rotation" (V x v ^ 0), and the 
corresponding inner caustic is a closed tube whose cross section is a section of the elliptic umbilic (£'-4) catastrophe, 
called a "caustic ring". The second case is characterized by irrotational flow (V X v = 0), which gives rise instead to 
inner caustics with a tent-like structure. Tidal torque theory predicts that ordinary cold dark matter such as WIMPs 
would be irrotational [10], hence giving rise to tentlike caustics. In contrast, "net overall rotation" (~ vortices) in an 
axionic BEC would result in ring caustics [5[ [7] . 

For the first case ("net overall rotation"), the radii of the caustic rings a n can be predicted in terms of a single 
parameter, j m ax, under the assumption that the time evolution of the galactic halo is self-similar llj. That is, the 
phase space structure of the halo is time-independent under appropriate rescaling of all distances, velocities, and 
densities [TlTU5| . If initially the overdensity in the halo has the profile 8M/M oc M~ e , where the parameter e is 
estimated to lie in the range 0.25 to 0.35 from the observed power spectrum of density perturbations on galaxy scales, 
the prediction for a n is given by [TTJ [TCI H7J 
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where v rot is the galactic rotation velocity and n = 1,2,3,... 

Intriguing hints for caustic rings with radii consistent with the prediction of Eq. have been found in a multitude 
of observations: the distribution of bumps in the rotation curve of the Milky Way (MW) [15], and the appearance 
of a triangular feature in the IRAS map of the MW in the precise direction tangent to the nearest caustic ring |18j . 
Outside our home galaxy, further corroboration can be found in the statistical distribution of bumps in a set of 32 
extended and well-measured galactic rotation curves [20] ; the recent rotation curve of the Andromeda galaxy, which 
shows three bumps at locations consistent with the prediction of 0; and gravitational lensing evidence in galaxy 
clusters [5T|. Together, these observations lend credence to the argument that an axion BEC constitutes a substantial 
fraction of dark matter. 

Further evidence that may favor the existence of caustic rings is the recently-discovered Monoceros Ring of stars, 
an over-density of stars in the plane of the Galaxy at a galactocentric distance of 20 kpc [251 H3] • Their proximity 
to the n — 2 inner caustic has led to speculations that their origin has to do with the presence of the ring [T5]. To 
test this claim, one will have to resort to numerical simulations involving the gravitational field near a caustic ring. 
Unfortunately, the expression for the gravitation field is in general in the form of an integral, whose evaluation is 
computationally expensive, rendering numerical simulations difficult. It is the goal of this paper to help circumvent 
this difficulty by evaluating the near-field caustic-ring gravitational field by analytical methods. 

This paper is structured as follows. In Section II, we present our calculation of the gravitational field near a caustic 
ring (in the limit of infinite radius), specializing to the case where £, a parameter which we shall define and expect to 
be close to one, equals unity. In Section III, we compare our analytic results with numerics, to establish the validity 
of our formulae in Section II. 

Finally, we note that, building on insights acquired in Section II, we can straightforwardly generalize our calculation 
to the case where £ 7^ 1. The computational method employed in Section II remains applicable, though the algebra 
becomes more cumbersome. Since the generalization is not particularly illuminating, it is omitted in this paper. 



II. CALCULATION 



A. Theoretical Framework 



Much of the discussion in this subsection is taken from |19j . Assuming axial symmetry and that the radius of the 
ring far exceeds the transverse dimensions of the caustic, the distribution of CDM near a caustic ring at some given 
time is given by 
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where p and z are the usual cylindrical coordinates [16j (the angular coordinates 4> is absent due to rotational 
symmetry). The parameters b,a,u,TQ, s characterize the caustic ring (see [19j for a good summary). The variable r 
is a temporal coordinate, and denotes the time when the CDM particle crosses the z — plane. Hence — r can be 
viewed as the age of the particle. The variable a = | — 8 is an angular coordinate, with 8 being the polar angle of the 
particle at its last turnaround. Hence we have two related sets of variables describing the dynamics of the particles: 
(p, z) and (a, r). The location of the caustic is obtained from Eqs. ^ and (3j) by solving for the values of a and r 
at which the Jacobian determinant between the two sets of coordinates vanishes: 
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The density of the dark matter particles in physical space is given by |19j 
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where otj(p,z) and Tj(p,z) are solutions to p(a,r) ~ p and z(a,r) = z. The number of solutions N(p,z) represents 



the number of flows at location (p,z). 
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- is the mass of dark matter particles falling in per unit solid angle and 
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For convenience, we introduce two new variables 



P = \ UT l ( 6 ) 



q = a nn. ' ( 7 ) 
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The variables p and g have clear physical meanings: p is the longitudinal dimension of the caustic (the caustic extends 
from atoa+pat z = 0), while g is the transverse dimension (the highest point of the caustic are at z — q/2). 
Under the assumption that p. q <C a (i.e. infinite caustic radius), the gravitational field g is given by |19j 

g(p,z) = -2G [ dp'dz'd(p',z') {P ~(^~ Z '\, 2 , (8) 

where p' is integrated from to oo and z' from — oo to oo. 

Performing a coordinate transformation from (p,z) to (a, r) and using Eq. ([5]), and ignoring the spatial variation 
of over the size of the caustic, we have 

2G dM f (p- p(a,r),z - z{a,r)) 

= / dadT l 7 v^T7 7 va- ( 9 ) 

p dildr J (p - p(a,T)y + (z - z(a,T)) z 

To simplify the calculation, we introduce the following rescaled variables 

T = — (10) 

A = * (11) 

UTq 

Z = p (12) 

X = (13) 
P 

and the dimensionless parameter £ = su/b 2 . The surface of the caustic is then given by: 

X s = (T-1)(2T-1) (14) 
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From Eqs. (14| and ( |T5| , we see that the caustic extends from X = —1/8 to X = 1, independent of the value of 

£. However, the highest (lowest) point of the caustic Z max ( Z mi „) does depend on £, and is given by ±|^/|. For 

£ = 1, the caustic possesses a Zj, symmetry, as its three cusps rotate into one another under a rotation by 120°. 
The gravitational field of the caustic ring has the form 

^-^mp^^W2S)) t (16) 

where I P {X, Z) and I Z {X, Z) are given respectively by 
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In the galactic plane (Z = 0) and assuming £ = 1, I z = by symmetry and the radial component I p has already 
been evaluated in IHfl: 
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In this paper, we aim to generalize this result to arbitrary X and Z. As in [115], we will restrict our analysis to 
the case £ = 1, which is of order its expected value. (For example, based on the IRAS map of the galactic disk, 
Sikivie estimated in [T5] that the n — 5 caustic ring has dimensions p ~ 130 pc and q ~ 200 pc. Using the relation 
p/q = -\/l6£/27, one then obtains £ ~ 0.71). The £ = 1 case has the additional advantage that the algebra is less 
complicated. Generalizing to the case £ ^ 1 is not difficult but rather tedious, and the method used in this paper 
applies straightforwardly. 



B. Gravitational field for { = 1 



Our goal in this section is to compute the following integrals 
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Eq. (j 1 9|) can be rewritten as 
dAdT 
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while Eq. (|20|) can be rewritten as 
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The poles c 3 , c 4 , C5, and C6 are simple, so the A integral can be evaluated straightforwardly using residue theory. 
The remaining task now is to determine the location of the poles at given X, Z, and T to see which contributes to 
the A integral. 

Using identities such as (c 3 — cs)(c 3 — cq) — 2(c 3 — ci)(c 3 — c 2 ) and (c 3 — C5)(c 3 — eg) = — |~(c3 — c 7)j we have 
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We generalize A to a complex variable, and use as our contour the real axis and the upper half-circle of infinite 
radius and counter-clockwise orientation. (Note that one could use instead the lower half-circle, and our procedure 
still works). The contribution from integrating over the upper-half-circle perimeter vanishes since our integrand has 
the form P(A)/Q(A), where P(A) and Q(A) are polynomials, and the degree of P(A) is less than that of Q(A) by two. 
Hence, each pole in the upper half plane would contribute 2iri times its residue to Eq. (19 1. Slight complication arises 



when the poles are on the real axis, in which case we will have to deform the contour to avoid the poles. However, 
the reality of the integral ensures that another pole enters the upper half plane simultaneously as another leaves it, 
so that our original recipe of each pole within the contour contributing 27ri x Residue remains valid. 

We now analyze the behavior of the poles as T changes, given fixed X and Z. The imaginary parts of (19 1 cj., C2, 
and C4 are 
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Im(c 3 ) = -T+\j\ (V^ 2 + (X-(T- 1) 2 ))2 + Z 2 + T 2 + (X - (T - 1)2)) = -Im(c 6 ) 
Im(c 4 ) = -T-\j\ (V^ 2 + (X-(T- 1)2))2 + Z 2 + T 2 + (X - (T - l) 2 ) = -Im(c 5 ) 
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We can solve for where (as a function of T) each pole crosses the real axis by setting any of the above to 0. The 
resultant condition is a quartic equation of the form 



4T 4 _ 8T 3 + 4 ^ _ X ^ T 2 _ Z 2 = 
To express the four roots of this equation in a more concise form, we first define the functions 
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With these, we further introduce 
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In the galactic plane Z = 0, Eq. p4| becomes T 4 - 2T 3 + (1 - A)T 2 = 0, and can be easily factored. The four 
roots are T = (twice), T = 1 ± -\/A. Hence, the gravitational held admits a particularly simple form in this case, as 
we alluded to earlier. 



When /4 vanishes, Eqs. (39 1 to (42 1 become undefined. This occurs when / 2 = and < 0. The first condition 
implies that Z = ±(1 - X)/\73, which means that the second condition becomes 

128(A - l) 3 + 1728Z 2 + 1152(A - f )Z 2 > 0, (43) 



which implies A > —1/8. Combining these conditions, we see that roots as given by Eqs. (39 1 to (42) are ill-defined 
on the two lines Z = ±(1 — A)/-\/3, X > —1/8. Note that the line corresponding to the minus sign connects the 
highest point (X — —1/8,Z — 3a/3/8) and the rightmost tip (X = 1, Z = 0) of the caustic, while the other is its 
vertical reflection about the Z = axis. 



To solve for the roots on these two lines, we use the facts A = and /3 < to first simplify Eq. (34), from which 
we can solve for the four roots: 

Ti = \ ( 1 - M +h ~ \i .A - h 2 X j ! /, ) j ( 44 ) 
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Ti = \ ( I 4- \/^f + ./>. - \/ .A - .A-2..V j ^ -/, j | (47) 



where 



f 5 (X) = ^{-l-6X + 15X 2 -8X 3 ) 1/3 . (48) 

Among the four roots, only the real ones correspond to poles crossing the real-axis of the A-complex plane. It turns 
out that qualitative differences exist between the behavior of the poles within and outside the caustic, so we have to 



study these two cases separately. Within the caustic (boundary given by (14)), all four roots are real; outside, only 
two of the four roots arc real. This difference has to do with why a kink in the gravitational field appears as we cross 
the boundary of the caustic ring, which we saw in the case of Z = discussed earlier. 

1. Behavior of the poles within the caustic 

Within the caustic, all four roots to the quartic equation are real, implying that the poles cross the real axis a total 
of four times. In particular, the pole c 3 would cross the real axis three times, while C4 would only do so once. Since 
only poles in the upper half plane contribute to the integral, we will have to partition our integral over T into multiple 
regions accordingly, and evaluate each partition separately. While it is nontrivial to determine precisely which of 
the roots of the quartic equation correspond to which pole, operationally this information is not essential, as we can 



always just reorder the roots and use them appropriately. That is, we compute the four roots using Eqs. (39) - (42) 
(Ti_i.4) and rank them in increasing order, and relabel them as T\ — > T4, where T\ < T% < T3 < T4. The smallest root 
T\ is the value of T at which C4 crosses the real axis (at fixed A and Z), while the other three larger roots are when 
C3 crosses the real axis. More specifically, for T < Ti, both C3 and C4 are in the upper half plane, and they would 
contribute 2iri x [Residue^) + Residue(c4)] to the integral. As T increases beyond Ti, C4 leaves the upper half and 
plane, and C5 enters, while C3 remains. Hence, now C3 and C5 would contribute 2ni x [Residue(c3) + Residue^)] to 
the integral. This holds until T = T2, when C3 also leaves the upper half plane and cq enters simultaneously, so now C5 
and C6 contribute to the integral. When T gets to T 3 , C3 reenters while c@ leaves. Finally, C3 leaves again for T > T 4 . 

Summarizing, at given A and Z, the poles that are in the upper half complex A plane as a function of T are given 
by: 
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• T e (-00, fi) : Im(c 3 ),Im(c 4 ) > 

• T G (Ti,f 2 ) : Im(c 3 ),Im(c 5 ) > 

• T e (f 2 ,f 3 ) : Im(c 5 ),Im(c 6 ) > 

• Te (T 3 ,T 4 ) : Im( C3 ),Im(c 5 ) > 

• T e (T 4 ,oo) : Im(c 5 ),Im(c 6 ) > 

Note also our calculation allows for an arbitrary constant which can be fixed by taking the limit that the caustic 
tube shrinks to a line. In this limit, the caustic is an infinite line at X = Z = 0, and the density of dark matter 
particles is proportional to 1/V 'X 2 + Z 2 . Using Gauss' law, we can show that I = f/2, where f is the radial unit 
vector pointing away from the origin (recall that I and g point in opposite directions). Hence, we have to add \ 
to our solution for P p and to P z . Alternatively, we could obtain these constants from the asymptotic behaviour of 
I p (£ = 1, X — > oo) at Z = (and taking advantage of the rotational symmetry that emerges in this limit of vanishing 
cross section size). 

Hence after integrating over A, the integral P p (inside caustic) becomes 

n - \ [jl +£) 1? {^ 2T+ l-X + ,Z + V-2T + ix-,z) + ? (49) 
which can be integrated to yield 



4P p = \y2T4, - 1 + X - iZ + \J2T 4 - 1 + X + iZj - 2T 3 - 1 + X - iZ + \J2T 3 - 1 + X + iZ 

(^2P 2 -\ + X -iZ + ^2f 2 -l + X + iZ^j - (^2^-1 + X -iZ + \J 2T i - 1 + X + iZ^j + 2 (50) 



This can be simplified further to 



2I l p = Re (^/2T 4 -1 + X -iZ - ^2T 3 - 1 + X - iZ + ^2T 2 - 1 + X - iZ - ^2T 1 - 1 + X - iZ + lj 



(51) 



Meanwhile, the integral P z becomes 



t 2 rT*\ / 1 1 



4 (X, + L 3 ) \^-2T +1- X + iZ V-2T+1-X -iZj ' 
which can be readily integrated to yield 



5 = 7 /. + L dT I (52) 



(-\fl - X - 2T 4 + iZ + \jl - X - 2T 4 - iZ^j + (y/l - X - 2T 3 + iZ - \j I - X - 2T 3 - i? 
-\Jl - X - 2T 2 + iZ + \Jl - X - 2T 2 - iZ^j + (^Jl - X - 2T X + iZ - \J I - X - 2T\ - iZ^j . (53) 
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As before, this can be simplified further to 

2P Z = Im ^yj 2T 4 — l + X — iZ— \j 2T 3 — l + X — iZ+\j 2T 2 — l + X — iZ— \j 2T\ — l + X — iZ + T^J. (54) 

Note that P p and P z are respectively the real and imaginary part of the same quantity. This is expected since the 
gravitational field is conservative. 



2. Behavior of the poles outside the caustic 



Outside the caustic tube, the quartic polynomial has two complex and two real solutions. Only the real solutions 
correspond to poles crossing the real axis. If we start with T very negative, only C3 and C4 are in the upper half plane 
(as Im(c5) = — Im(c 3 ) and Im(c6) = — Im(c 4 )). As T increases, first c 4 crosses the real axis and leaves the upper 
half complex plane while C5 enters. After that C3 also exits while cq enters simultaneously. We denote the two real 



solutions to Eqs. (39) to (|42|) as T\ and T 2 , with T± < T 2 . At given X and Z, the poles that are in the upper half 



complex plane are given by 

• T G (-00, fx) : Im(c 3 ),Im(c 4 ) > 

• TG (fi.fa) :Im(c3),Im(c 5 ) >0 

• TG (T 2 ,oo) : Im(c 5 ),Im(c 6 ) > 0. 



Using residue theory, one obtains as before that the radial component of I outside the caustic (as before, we add ~ 



to 1°) 



41° 

p 



2T 2 -1 + X -iZ+\l 2T 2 - 1 + x + ;z 
Again, one can simplify this to 



The z component of I is 



21° = Re ( \f2T 2 - 1 + X - iZ - \f2Ti - 1 + X - iZ + 1 



2Ti - 1 + V2Ti - 1 + X + +2. (55) 



(56) 



41° = -y 2T 2 - 1 + X - iZ + V2T 2 - 1 + X + iZ + \/2Tx - 1 + X - iZ - \f2Tx - 1 + X + iZ, 
which can be simplified to 



■21.' = J ni j \, 2T 2 - 1 + X-iZ- \j2Ti -1 + X-iZ + l 
We plot the vector field / calculated using Eqs. (51), (54), (56), (58) in Figure [l] 



(57) 



(58) 



III. COMPARISON WITH NUMERICS 



To check the validity of our computation, we compare results obtained by our formulae with numerics. For better 
convergence in our numerical calculation, we rewrite Eqs. (19) and (20) in polar coordinates (A = Rcos(6); T = 
l + Rsm(6)), 



dRdO 



Rcos{26)[X - R 2 cos(20)] 



2ir [X - R 2 cos(2<2)] 2 + [Z- 2 J Rcos(0)(l + Rsin(6))} 2 ' 
dRd6 Rcos(26)[Z - 2i?cos(6>)(l + i?sin(6»))] 

2tt [X - R 2 cos(26»)] 2 + [Z- 2Rcos(6)(l + R sin(6»))] 2 ' 



(59) 
(60) 



and first perform the integration in the angular coordinate 9 in our numerical routine. The R integral is evaluated 
from to 100. 

The analytical calculation is carried out using a C- 



computer program, based on Eqs. (51), (54), (56), (58). Due 



to the presence of floating-point errors, we have to introduce a threshold value in our program to decide whether a 
root to Eq.(34) is real or not: a root is considered real if its imaginary part is less than 1 x 10 -5 . 



The comparison between analytics and numerics has been checked extensively at different values of X and Z, and 
is found to be in excellent agreement in every case, as is shown in Fig. [2j for Z = 0.1, 0.5, 1 and — 1 < X < 2. 
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FIG. 1: Vector representation of —I(X, Z) (not g) near a caustic ring. A silhouette of the caustic cross section is 
plotted in the background by plotting parametrically Eqs. ^ and ([3]). The curves in the background are contours 

of constant a and r. 



IV. CONCLUSIONS 



The presence of cold and collisionless dark matter particles entails the formation of dark matter caustics in galactic 
halos. Evidence for inner caustic rings (as opposed to tents) exist, rendering support for the case that dark matter 
particles possess angular momentum and undergo net overall rotation. This can be explained if Bose-condensed axions 
constitute a substantial portion of dark matter. Due to their gravitational effects, caustic rings inevitably play an 
important role in the formation of structure in the early universe. However, incorporating them into numerical simu- 
lations of structure formation is nontrivial, because the gravitational field of a caustic ring is traditionally expressed 
in the form of an integral, whose evaluation is computationally very expensive. As a first step towards overcoming 
this difficulty, we evaluate in this paper the integral analytically in the vicinity of the ring (equivalently, taking the 
limit of infinite radius) . Our calculation is based primarily on the use of residue theory, in which non-analyticity in 
the gravitational field corresponds to poles crossing the real axis. We find that our analytical formulae are in good 
agreement with numerics. An interesting subsequent project would be to generalize our calculation to caustic rings 
of finite radius. 
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FIG. 2: I r and I z for — 1 < X < 2 at Z — 0.1, 0.5, 1. The continuous curve is obtained by evaluating the analytical 
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